Signal Decay¶
In this chapter, we use simulated data to show how BOLD signal decays and how we can glean useful information from that decay rate.
import os
import imageio
import matplotlib.pyplot as plt
import nibabel as nib
import numpy as np
import seaborn as sns
from IPython import display
from matplotlib.animation import FuncAnimation
from myst_nb import glue
from nilearn import image, plotting
from nilearn.glm import first_level
from repo2data.repo2data import Repo2Data
from scipy import signal
from book_utils import predict_bold_signal
sns.set_style("whitegrid")
plt.rcParams.update({
"text.usetex": True,
"font.family": "sans-serif",
"font.sans-serif": ["Helvetica"]
})
# Install the data if running locally, or point to cached data if running on neurolibre
DATA_REQ_FILE = os.path.join("../binder/data_requirement.json")
# Download data
repo2data = Repo2Data(DATA_REQ_FILE)
data_path = repo2data.install()
data_path = os.path.abspath(os.path.join(data_path[0], "data"))
out_dir = os.path.join(data_path, "signal-decay")
os.makedirs(out_dir, exist_ok=True)
---- repo2data starting ----
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/repo2data
Config from file :
../binder/data_requirement.json
Destination:
./../data/multi-echo-data-analysis
Info : ./../data/multi-echo-data-analysis already downloaded
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/nilearn/datasets/__init__.py:96: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
"Numpy arrays.", FutureWarning)
Signal decays as echo time increases¶
func_dir = os.path.join(data_path, "sub-04570/func/")
data_files = [
os.path.join(func_dir, "sub-04570_task-rest_echo-1_space-scanner_desc-partialPreproc_bold.nii.gz"),
os.path.join(func_dir, "sub-04570_task-rest_echo-2_space-scanner_desc-partialPreproc_bold.nii.gz"),
os.path.join(func_dir, "sub-04570_task-rest_echo-3_space-scanner_desc-partialPreproc_bold.nii.gz"),
os.path.join(func_dir, "sub-04570_task-rest_echo-4_space-scanner_desc-partialPreproc_bold.nii.gz"),
]
echo_times = [12., 28., 44., 60.]
img = image.index_img(data_files[0], 0)
data = img.get_fdata()
vmax = np.max(data)
idx = np.vstack(np.where(data > 1000))
min_x, min_y, min_z = np.min(idx, axis=1)
max_x, max_y, max_z = np.max(idx, axis=1)
imgs = []
for f in data_files:
img = image.index_img(f, 0)
data = img.get_fdata()
data = data[min_x:max_x, min_y:max_y, min_z:max_z]
img = nib.Nifti1Image(data, img.affine, img.header)
imgs.append(img)
plt.style.use("dark_background")
fig, axes = plt.subplots(figsize=(26, 4), ncols=len(data_files))
for i_echo, img in enumerate(imgs):
te = echo_times[i_echo]
if i_echo == 0:
data = img.get_fdata()
vmax = np.max(data)
plotting.plot_epi(
img,
cut_coords=[0],
display_mode="z",
annotate=False,
draw_cross=False,
axes=axes[i_echo],
figure=fig,
vmax=vmax,
cmap="gray",
)
axes[i_echo].set_title(f"TE={te}ms", fontsize=20, pad=0)
glue("fig_brain_decay", fig, display=False)
# Reset the style
plt.style.use("default")
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
253 command, cwd=cwd if cwd is not None else self.texcache,
--> 254 stderr=subprocess.STDOUT)
255 except FileNotFoundError as exc:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in check_output(timeout, *popenargs, **kwargs)
410 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411 **kwargs).stdout
412
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
487
--> 488 with Popen(*popenargs, **kwargs) as process:
489 try:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)
799 errread, errwrite,
--> 800 restore_signals, start_new_session)
801 except:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
1550 err_msg += ': ' + repr(err_filename)
-> 1551 raise child_exception_type(errno_num, err_msg, err_filename)
1552 raise child_exception_type(err_msg)
FileNotFoundError: [Errno 2] No such file or directory: 'latex': 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
339 pass
340 else:
--> 341 return printer(obj)
342 # Finally look for special method names
343 method = get_real_method(obj, self.print_method)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
149 FigureCanvasBase(fig)
150
--> 151 fig.canvas.print_figure(bytes_io, **kw)
152 data = bytes_io.getvalue()
153 if fmt == 'svg':
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2228 else suppress())
2229 with ctx:
-> 2230 self.figure.draw(renderer)
2231
2232 if bbox_inches:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
72 @wraps(draw)
73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74 result = draw(artist, renderer, *args, **kwargs)
75 if renderer._rasterizing:
76 renderer.stop_rasterizing()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
2779 self.patch.draw(renderer)
2780 mimage._draw_list_compositing_images(
-> 2781 renderer, self, artists, self.suppressComposite)
2782
2783 for sfig in self.subfigs:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/_api/deprecation.py in wrapper(*inner_args, **inner_kwargs)
429 else deprecation_addendum,
430 **kwargs)
--> 431 return func(*inner_args, **inner_kwargs)
432
433 return wrapper
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
2879 artists.remove(spine)
2880
-> 2881 self._update_title_position(renderer)
2882
2883 if not self.axison or inframe:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in _update_title_position(self, renderer)
2810 if (ax.xaxis.get_ticks_position() in ['top', 'unknown']
2811 or ax.xaxis.get_label_position() == 'top'):
-> 2812 bb = ax.xaxis.get_tightbbox(renderer)
2813 else:
2814 bb = ax.get_window_extent(renderer)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in get_tightbbox(self, renderer, for_layout_only)
1081 ticks_to_draw = self._update_ticks()
1082
-> 1083 self._update_label_position(renderer)
1084
1085 # go back to just this axis's tick labels
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _update_label_position(self, renderer)
2078 # get bounding boxes for this axis and any siblings
2079 # that have been set by `fig.align_xlabels()`
-> 2080 bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
2081
2082 x, y = self.label.get_position()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_boxes_siblings(self, renderer)
1866 axis = ax._get_axis_map()[axis_name]
1867 ticks_to_draw = axis._update_ticks()
-> 1868 tlb, tlb2 = axis._get_tick_bboxes(ticks_to_draw, renderer)
1869 bboxes.extend(tlb)
1870 bboxes2.extend(tlb2)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_bboxes(self, ticks, renderer)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in <listcomp>(.0)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in get_window_extent(self, renderer, dpi)
901
902 with cbook._setattr_cm(self.figure, dpi=dpi):
--> 903 bbox, info, descent = self._get_layout(self._renderer)
904 x, y = self.get_unitless_position()
905 x, y = self.get_transform().transform((x, y))
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in _get_layout(self, renderer)
306 _, lp_h, lp_d = renderer.get_text_width_height_descent(
307 "lp", self._fontproperties,
--> 308 ismath="TeX" if self.get_usetex() else False)
309 min_dy = (lp_h - lp_d) * self._linespacing
310
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
228 fontsize = prop.get_size_in_points()
229 w, h, d = texmanager.get_text_width_height_descent(
--> 230 s, fontsize, renderer=self)
231 return w, h, d
232
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
397 else:
398 # use dviread.
--> 399 dvifile = self.make_dvi(tex, fontsize)
400 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
401 page, = dvi
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in make_dvi(self, tex, fontsize)
291 self._run_checked_subprocess(
292 ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 293 texfile], tex, cwd=tmpdir)
294 (Path(tmpdir) / Path(dvifile).name).replace(dvifile)
295 return dvifile
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
256 raise RuntimeError(
257 'Failed to process string with tex because {} could not be '
--> 258 'found'.format(command[0])) from exc
259 except subprocess.CalledProcessError as exc:
260 raise RuntimeError(
RuntimeError: Failed to process string with tex because latex could not be found
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
253 command, cwd=cwd if cwd is not None else self.texcache,
--> 254 stderr=subprocess.STDOUT)
255 except FileNotFoundError as exc:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in check_output(timeout, *popenargs, **kwargs)
410 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411 **kwargs).stdout
412
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
487
--> 488 with Popen(*popenargs, **kwargs) as process:
489 try:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)
799 errread, errwrite,
--> 800 restore_signals, start_new_session)
801 except:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
1550 err_msg += ': ' + repr(err_filename)
-> 1551 raise child_exception_type(errno_num, err_msg, err_filename)
1552 raise child_exception_type(err_msg)
FileNotFoundError: [Errno 2] No such file or directory: 'latex': 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
339 pass
340 else:
--> 341 return printer(obj)
342 # Finally look for special method names
343 method = get_real_method(obj, self.print_method)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
149 FigureCanvasBase(fig)
150
--> 151 fig.canvas.print_figure(bytes_io, **kw)
152 data = bytes_io.getvalue()
153 if fmt == 'svg':
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2228 else suppress())
2229 with ctx:
-> 2230 self.figure.draw(renderer)
2231
2232 if bbox_inches:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
72 @wraps(draw)
73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74 result = draw(artist, renderer, *args, **kwargs)
75 if renderer._rasterizing:
76 renderer.stop_rasterizing()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
2779 self.patch.draw(renderer)
2780 mimage._draw_list_compositing_images(
-> 2781 renderer, self, artists, self.suppressComposite)
2782
2783 for sfig in self.subfigs:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
49 renderer.start_filter()
50
---> 51 return draw(artist, renderer, *args, **kwargs)
52 finally:
53 if artist.get_agg_filter() is not None:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/_api/deprecation.py in wrapper(*inner_args, **inner_kwargs)
429 else deprecation_addendum,
430 **kwargs)
--> 431 return func(*inner_args, **inner_kwargs)
432
433 return wrapper
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
2879 artists.remove(spine)
2880
-> 2881 self._update_title_position(renderer)
2882
2883 if not self.axison or inframe:
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in _update_title_position(self, renderer)
2810 if (ax.xaxis.get_ticks_position() in ['top', 'unknown']
2811 or ax.xaxis.get_label_position() == 'top'):
-> 2812 bb = ax.xaxis.get_tightbbox(renderer)
2813 else:
2814 bb = ax.get_window_extent(renderer)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in get_tightbbox(self, renderer, for_layout_only)
1081 ticks_to_draw = self._update_ticks()
1082
-> 1083 self._update_label_position(renderer)
1084
1085 # go back to just this axis's tick labels
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _update_label_position(self, renderer)
2078 # get bounding boxes for this axis and any siblings
2079 # that have been set by `fig.align_xlabels()`
-> 2080 bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
2081
2082 x, y = self.label.get_position()
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_boxes_siblings(self, renderer)
1866 axis = ax._get_axis_map()[axis_name]
1867 ticks_to_draw = axis._update_ticks()
-> 1868 tlb, tlb2 = axis._get_tick_bboxes(ticks_to_draw, renderer)
1869 bboxes.extend(tlb)
1870 bboxes2.extend(tlb2)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_bboxes(self, ticks, renderer)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in <listcomp>(.0)
1062 """Return lists of bboxes for ticks' label1's and label2's."""
1063 return ([tick.label1.get_window_extent(renderer)
-> 1064 for tick in ticks if tick.label1.get_visible()],
1065 [tick.label2.get_window_extent(renderer)
1066 for tick in ticks if tick.label2.get_visible()])
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in get_window_extent(self, renderer, dpi)
901
902 with cbook._setattr_cm(self.figure, dpi=dpi):
--> 903 bbox, info, descent = self._get_layout(self._renderer)
904 x, y = self.get_unitless_position()
905 x, y = self.get_transform().transform((x, y))
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in _get_layout(self, renderer)
306 _, lp_h, lp_d = renderer.get_text_width_height_descent(
307 "lp", self._fontproperties,
--> 308 ismath="TeX" if self.get_usetex() else False)
309 min_dy = (lp_h - lp_d) * self._linespacing
310
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
228 fontsize = prop.get_size_in_points()
229 w, h, d = texmanager.get_text_width_height_descent(
--> 230 s, fontsize, renderer=self)
231 return w, h, d
232
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
397 else:
398 # use dviread.
--> 399 dvifile = self.make_dvi(tex, fontsize)
400 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
401 page, = dvi
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in make_dvi(self, tex, fontsize)
291 self._run_checked_subprocess(
292 ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 293 texfile], tex, cwd=tmpdir)
294 (Path(tmpdir) / Path(dvifile).name).replace(dvifile)
295 return dvifile
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
256 raise RuntimeError(
257 'Failed to process string with tex because {} could not be '
--> 258 'found'.format(command[0])) from exc
259 except subprocess.CalledProcessError as exc:
260 raise RuntimeError(
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 1872x288 with 8 Axes>
<Figure size 1872x288 with 8 Axes>
Fig. 1 Signal decay in the brain.¶
The impact of \(T_{2}^{*}\) and \(S_{0}\) fluctuations on BOLD signal¶
# Simulate data
MULTIECHO_TES = np.array([15, 30, 45, 60, 75, 90])
SINGLEECHO_TE = np.array([30])
# For a nice, smooth curve
FULLCURVE_TES = np.arange(0, 101, 1)
n_echoes = len(FULLCURVE_TES)
pal = sns.color_palette('cubehelix', 8)
Simulate \(T_{2}^{*}\) and \(S_{0}\) fluctuations¶
# Simulate data
# We'll convolve with HRF just for smoothness
hrf = first_level.spm_hrf(1, oversampling=1)
N_VOLS = 21
SCALING_FRACTION = 0.1 # used to scale standard deviation
MEAN_T2S = 35
MEAN_S0 = 16000
# simulate the T2*/S0 time series
# The original time series will be a random time series from a normal distribution,
# convolved with the HRF
ts = np.random.normal(loc=0, scale=1, size=(N_VOLS+20,))
ts = signal.convolve(ts, hrf)[20:N_VOLS+20]
ts *= SCALING_FRACTION / np.std(ts)
ts -= np.mean(ts)
t2s_ts = (ts * MEAN_T2S) + MEAN_T2S
s0_ts = (ts * MEAN_S0) + MEAN_S0
Plot BOLD signal decay for a standard single-echo scan¶
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [30])
fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
)
ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
glue("fig_signal_decay_single-echo", fig, display=False)
Fig. 2 BOLD signal decay for a standard single-echo scan¶
Plot BOLD signal decay and BOLD contrast¶
fullcurve_signal_active = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [40])
fullcurve_signal_inactive = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [20])
fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
FULLCURVE_TES,
fullcurve_signal_active[:, 0],
alpha=0.5,
color="red",
label="High Activity",
)
ax.plot(
FULLCURVE_TES,
fullcurve_signal_inactive[:, 0],
alpha=0.5,
color="blue",
label="Low Activity",
)
ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.legend(fontsize=20)
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
glue("fig_signal_decay_contrast", fig, display=False)
Fig. 3 BOLD signal decay and BOLD contrast¶
Plot single-echo data resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how fMRI data fluctuates over time.
Plot single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how single-echo data is a sample from a signal decay curve.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
# gif code based on https://www.geeksforgeeks.org/create-an-animated-gif-using-python-matplotlib/
fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})
# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
ax0_line_plot = axes[0].plot(
singleecho_signal[0, :],
color="black",
zorder=0
)
ax0_scatter_plot = axes[0].scatter(
0,
singleecho_signal[:, 0],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_line_plot = axes[1].plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
zorder=0,
)[0]
ax1_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, 0],
color="orange",
s=150,
alpha=1.,
label="Single-Echo Signal",
zorder=1,
)
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack((
frame,
singleecho_signal[:, frame],
))
)
ax1_line_plot.set_data((
FULLCURVE_TES,
fullcurve_signal[:, frame],
))
ax1_scatter_plot.set_offsets(
np.column_stack((
SINGLEECHO_TE,
singleecho_signal[:, frame],
))
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay", html, display=False)
Fig. 4 Single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how changes in fMRI data can be driven by both S0 and T2* fluctuations.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})
t2s_value = predict_bold_signal(
np.array([t2s_ts[0]]),
np.array([s0_ts[0]]),
np.array([t2s_ts[0]])
)[0]
ax0_line_plot = axes[0].plot(
singleecho_signal[0, :],
color="black",
zorder=0
)
ax0_scatter_plot = axes[0].scatter(
0,
singleecho_signal[:, 0],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_line_plot = axes[1].plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
zorder=0,
)[0]
ax1_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, 0],
color="orange",
s=150,
alpha=1.,
label="Single-Echo Signal",
zorder=1,
)
ax1_t2s_scatter_plot = axes[1].scatter(
t2s_ts[0],
t2s_value,
marker="*",
color="blue",
s=500,
alpha=0.5,
label="$T_{2}^{*}$",
zorder=1,
)
ax1_s0_scatter_plot = axes[1].scatter(
0,
s0_ts[0],
marker="*",
color="red",
s=500,
alpha=0.5,
label="$S_{0}$",
clip_on=False,
zorder=2,
)
# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
t2s_value = predict_bold_signal(
np.array([t2s_ts[frame]]),
np.array([s0_ts[frame]]),
np.array([t2s_ts[frame]])
)[0]
ax0_scatter_plot.set_offsets(
np.column_stack((
frame,
singleecho_signal[:, frame],
))
)
ax1_line_plot.set_data((
FULLCURVE_TES,
fullcurve_signal[:, frame],
))
ax1_scatter_plot.set_offsets(
np.column_stack((
SINGLEECHO_TE,
singleecho_signal[:, frame],
))
)
ax1_t2s_scatter_plot.set_offsets(
np.column_stack((
t2s_ts[frame],
t2s_value,
))
)
ax1_s0_scatter_plot.set_offsets(
np.column_stack((
0,
s0_ts[frame],
))
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay2", html, display=False)
Fig. 5 Single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how fluctuations in S0 and T2* produce different patterns in the full signal decay curves.
s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})
t2s_value = predict_bold_signal(
np.array([t2s_ts[0]]),
np.array([s0_ts[0]]),
np.array([t2s_ts[0]])
)[0]
ax0_line_plot = axes[0].plot(
ts,
color="black",
zorder=0
)
ax0_scatter_plot = axes[0].scatter(
0,
ts[0],
color="purple",
s=150,
label="Single-Echo Signal",
zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, 0],
alpha=0.5,
color="red",
label="$S_0$-Driven Decay Curve",
)[0]
ax1_t2s_line_plot = axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, 0],
alpha=0.5,
color="blue",
label="$T_{2}^{*}$-Driven Decay Curve",
)[0]
# Set constant params for figure
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack((
frame,
ts[frame],
))
)
ax1_t2s_line_plot.set_data((
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, frame],
))
ax1_s0_line_plot.set_data((
FULLCURVE_TES,
s0based_fullcurve_signal[:, frame],
))
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay3", html, display=False)
Fig. 6 \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data¶
This shows how single-echo data, on its own, cannot distinguish between S0 and T2* fluctuations.
s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
s0based_singleecho_signal = s0based_fullcurve_signal[SINGLEECHO_TE, :]
t2sbased_singleecho_signal = t2sbased_fullcurve_signal[SINGLEECHO_TE, :]
fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})
ax0_line_plot = axes[0].plot(
ts,
color="black",
zorder=0
)
ax0_scatter_plot = axes[0].scatter(
0,
ts[0],
color="purple",
s=150,
zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, 0],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
zorder=0,
)[0]
ax1_t2s_line_plot = axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, 0],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
zorder=1,
)[0]
ax1_s0_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
s0based_singleecho_signal[:, 0],
color="red",
s=150,
alpha=1.,
label="$S_0$-Driven Single-Echo Signal",
zorder=2,
)
ax1_t2s_scatter_plot = axes[1].scatter(
SINGLEECHO_TE,
t2sbased_singleecho_signal[:, 0],
color="blue",
s=150,
alpha=1.,
label="$T_{2}^{*}$-Driven Single-Echo Signal",
zorder=3,
)
# Set constant params for figure
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack((
frame,
ts[frame],
))
)
ax1_s0_line_plot.set_data((
FULLCURVE_TES,
s0based_fullcurve_signal[:, frame],
))
ax1_t2s_line_plot.set_data((
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, frame],
))
ax1_s0_scatter_plot.set_offsets(
np.column_stack((
SINGLEECHO_TE,
s0based_singleecho_signal[:, frame],
))
)
ax1_t2s_scatter_plot.set_offsets(
np.column_stack((
SINGLEECHO_TE,
t2sbased_singleecho_signal[:, frame],
))
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay4", html, display=False)
Fig. 7 \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data¶
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data¶
This shows how S0 and T2* fluctuations produce different patterns in multi-echo data.
s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
s0based_multiecho_signal = s0based_fullcurve_signal[MULTIECHO_TES, :]
t2sbased_multiecho_signal = t2sbased_fullcurve_signal[MULTIECHO_TES, :]
fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})
ax0_line_plot = axes[0].plot(
ts,
color="black",
zorder=0,
)[0]
ax0_scatter_plot = axes[0].scatter(
0,
ts[0],
color="purple",
s=150,
zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, 0],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
zorder=0,
)[0]
ax1_t2s_line_plot = axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, 0],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
zorder=1,
)[0]
ax1_s0_scatter_plot = axes[1].scatter(
MULTIECHO_TES,
s0based_multiecho_signal[:, 0],
color="red",
s=150,
alpha=1.,
label="$S_0$-Driven Multi-Echo Signal",
zorder=2,
)
ax1_t2s_scatter_plot = axes[1].scatter(
MULTIECHO_TES,
t2sbased_multiecho_signal[:, 0],
color="blue",
s=150,
alpha=1.,
label="$T_{2}^{*}$-Driven Multi-Echo Signal",
zorder=3,
)
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
def AnimationFunction(frame):
"""Function takes frame as an input."""
ax0_scatter_plot.set_offsets(
np.column_stack((
frame,
ts[frame],
))
)
ax1_s0_line_plot.set_data((
FULLCURVE_TES,
s0based_fullcurve_signal[:, frame],
))
ax1_t2s_line_plot.set_data((
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, frame],
))
ax1_s0_scatter_plot.set_offsets(
np.column_stack((
MULTIECHO_TES,
s0based_multiecho_signal[:, frame],
))
)
ax1_t2s_scatter_plot.set_offsets(
np.column_stack((
MULTIECHO_TES,
t2sbased_multiecho_signal[:, frame],
))
)
anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay5", html, display=False)
Fig. 8 \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data¶
Plot \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)¶
When the BOLD signal is driven entirely by T2* fluctuations, the BOLD signal is very similar to a scaled version of those T2* fluctuations, but there are small differences.
fig, ax = plt.subplots(figsize=(16, 4))
scalar = np.linalg.lstsq(
t2s_ts[:, None],
t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
rcond=None
)[0]
ax.plot(
t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
color="orange",
label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
t2s_ts * scalar,
label="$T_{2}^{*}$ (scaled to signal)",
linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
glue("fig_t2s_bold_single-echo", fig, display=False)
Fig. 9 \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)¶
Plot \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)¶
When the BOLD signal is driven entirely by S0 fluctuations, the BOLD signal ends up being a scaled version of the S0 fluctuations.
fig, ax = plt.subplots(figsize=(16, 4))
scalar = np.linalg.lstsq(
s0_ts[:, None],
s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
rcond=None
)[0]
ax.plot(
s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
color="orange",
label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
s0_ts * scalar,
label="S0 (scaled to signal)",
linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
glue("fig_s0_bold_single-echo", fig, display=False)
Fig. 10 \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)¶